#include <iostream>
#include <fstream>
#include <stdlib.h>
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TFile.h"
#include <cmath>
#include "TGraph.h"
#include "TH2F.h"
#include "TTree.h"
#include <sstream>

using namespace std;

int main(int argc, char* argv[]) {
	//TApplication theApp("theApp", 0, argv); //line to allow interactive canvas in the program

	string ifname(argv[1]);
	string ofname = argv[2];
	string treename = argv[3];
	double d1, d2;
	d1 = d2 = 0;
	d1 = atof(argv[4]);
	d2 = atof(argv[5]);

	cout << "reading data from " << ifname << endl;
	ifstream ifile(ifname.c_str());
	if (!ifile) {
		cout << "file " << ifname << " not found" << endl;
		exit(1);
	};
	TFile *ofile = new TFile(ofname.c_str(), "RECREATE");
	const int nVariables = 24;
	string linenames[nVariables] = { "K", "GEN", "Weight", "T", "X",
				"Y", "S", "E", "Px", "Py", "Ps", "Sx",
				"Sy", "Ss", "theta", "phi", "x1", "y1", "z1", "r1", "x2", "y2",
				"z2", "r2" };
//	string linenames[nVariables] = { "K", "GEN", "Weight", "T(m)", "X(m)",
//						"Y(m)", "S(m)", " E(eV)", "Px(eV_c)", "Py(eV_c)", "Ps(eV_c)", "Sx",
//						"Sy", "Ss", "theta", "phi", "x1", "y1", "z1", "r1", "x2", "y2",
//						"z2", "r2" };
	Float_t data[nVariables];

	TTree *tree = new TTree(treename.c_str(), "data from file");

	for (int i = 0; i < nVariables; i++) {
		stringstream ssname2;
		ssname2 << linenames[i] << "/F";
		tree->Branch(linenames[i].c_str(), &data[i], ssname2.str().c_str());
	};
	double counter;
	counter = 0;
	double theta;
	double phi;
	double r1, r2, x1, x2, y1, y2, z1, z2;
	theta = phi = r1 = r2 = x1 = x2 = y1 = y2 = z1 = z2 = 0;
	double x0, y0, z0, e, px0, py0, pz0;
	x0 = y0 = z0 = e = px0 = py0 = pz0 = 0;
	//while (ifile&&counter<1000) {
		while (ifile) {
				counter++;
		for (int i = 0; i < 14; i++) {
			ifile >> data[i];
			if (i == 7)
				data[i] = data[i] / 1000.; //energy in keV for graphs
		};
		x0 = data[4];
		y0 = data[5];
		z0 = data[6];
		e = data[7];
		px0 = data[8];
		py0 = data[9];
		pz0 = data[10];

		theta = sqrt(px0 * px0 + py0 * py0);
		theta = atan2(theta, pz0);
		data[14] = theta;
		data[15] = phi = atan2(py0, px0);
		//double t1=
		data[16] = x1 = x0 + d1 * px0/pz0;
		data[17] = y1 = y0 + d1 * py0/pz0;
		data[18] = z1 = d1;
		data[19] = r1 = sqrt(x1 * x1 + y1 * y1);
		data[20] = x2 = x0 + d2 * px0/pz0;
		data[21] = y2 = y0 + d2 * py0/pz0;
		data[22] = z2 = d2;
		data[23] = r2 = sqrt(x2 * x2 + y2 * y2);
		//for(int i=0;i<nVariables;i++) cout<<data[i]<<" ";
		//cout<<endl;
		tree->Fill();
	};
	tree->Write("", TObject::kOverwrite);
	ofile->Close();
}
;
